﻿#include "IsoVortex.h"
#include "Mesh.h"
#include "worklet/ComputeResidual.h"
#include "worklet/Context.h"
#include "worklet/InitialCondition.h"
#include "worklet/MeshModifier.h"

struct IsoVortex::Context : EquationContext
{
  Context(const Configuration& cfg)
    : EquationContext(cfg)
  {
  }


  VTKM_EXEC ValueType ExactSolution(const FloatDefault t, const Vec3f& p) const
  {
    ValueType value;
    Vec3f _vel{ 1, 1, 0 };
    FloatDefault x = p[0] - _vel[0] * t;
    FloatDefault y = p[1] - _vel[1] * t;
    FloatDefault z = p[2] - _vel[2] * t;
    FloatDefault gam = 1.4, gamm1 = gam - 1, epi = 5.;
    FloatDefault xb, yb, r2;

    FloatDefault PI = vtkm::Pi();
    xb = x + 5;
    yb = y + 5;
    r2 = xb * xb + yb * yb;

    FloatDefault u = _vel[0] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * (-yb);
    FloatDefault v = _vel[1] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * xb;
    FloatDefault T = 1. - gamm1 * epi * epi / (8 * gam * PI * PI) * exp(1. - r2);
    FloatDefault density = pow(T, 1 / gamm1);
    FloatDefault pre = pow(density, gam);

    value[0] = density;
    value[1] = density * u;
    value[2] = density * v;
    value[3] = density * 0;
    value[4] = pre / gamm1 + 0.5 * density * (u * u + v * v);

    return value;
  }

  template<typename CoordType>
  VTKM_EXEC_CONT ValueType InitialCondition(const CoordType& p)
  {
    auto exact = ExactSolution(0.0, p);
    return exact;
    //return PrimaryToConservation(exact);
  }

  template<typename CoordType>
  VTKM_EXEC_CONT ValueType
  BoundaryCondition(const BCType& bc_type__, const ValueType& ul, const Vec3f& normal, const CoordType& p)
  {
    ValueType ur;
    BCEnum bc_type = static_cast<BCEnum>(bc_type__);
    ur = ExactSolution(0.0, p);
    return ur;
  }

  Context PrepareForExecution(vtkm::cont::DeviceAdapterId, vtkm::cont::Token&) const
  {
    return *this;
  }
};

IsoVortex::IsoVortex(const Configuration& cfg)
  : EulerEquation(cfg)
  , _context(new Context(cfg))
  , _finite_volume(Get<bool>("fv", false))
{
}

void IsoVortex::InitialCondition()
{
  _timer.Start();
  //Loops::InitialConditionIsoVortex(*this);
  std::cout << _mesh.GetCoordinate().GetNumberOfPoints() << std::endl;

  vtkm::cont::Invoker{}(Loops::WorkletInitCondition{},
                        _mesh.GetCellSet(),
                        _solution,
                        _mesh.GetCoordinate().GetDataAsMultiplexer(),
                        *_context);

  console::Success("IsoVortex::InitialCondition初始条件设置完成", _timer.GetElapsedTime(), "s");
}

void IsoVortex::UpdateResidual()
{
  _timer.Start();
  if (_finite_volume)
  {
    auto& transient = static_cast<Transient&>(*(_app.GetExecutioner()));
    Real current_time = transient.CurrentTime();

    vtkm::cont::Invoker{}(Loops::FiniteVolumeWorklet{ current_time },
                          _mesh.GetCellSet(),
                          _mesh.GetCoordinate(),
                          _mesh.GetCellLinks(),
                          this->_solution,
                          this->_rhs,
                          *this->_context);

    console::Info("IsoVortex::有限体积，", "时间: ", _timer.GetElapsedTime(), " s");
  }
  //else
  //{
  //  Loops::ComputeResidualIsoVortex(*this);
  //  console::Info("IsoVortex::有限差分，", "时间: ", _timer.GetElapsedTime(), " s");
  //}
}
